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1. 


INTRODUCTION 


This  report  documents  the  work  performed  during  a  Phase  I  SBIR  study  sponsored 
by  AFOSR.  The  next  phase  of  work  (Phase  H)  is  expected  to  be  two  years  in  duration. 
The  objective  of  the  overall  project  (Phases  I  and  II)  is  to  develop  an 
experimental /computational  approach  to  model  precombustion  chemistry  and 
transport  behavior  of  hydrocarbon  fuels  under  supercritical  conditions. 

The  Phase  I  study  focused  on  addressing  heat  transfer  characteristics  under  critical 
and  supercritical  conditions.  The  model  development  and  computations  were 
performed  at  CFDRC  by  the  principal  investigator  and  parallel  experimental  work 
was  conducted  by  Professor  L.  D.  Chen  at  the  University  of  Iowa  (under  a  subcontract 
from  CFDRC). 

1.1  Need  for  Predicting  Supercritical  Behavior  of  Fuels 

High  performance,  high  speed  aircraft  of  the  future  will  require  onboard  fuel  to  act 
as  a  heat  sink  to  meet  the  thermal  loads  arising  from  sources  such  as  the  engine,  the 
avionics,  the  air  frame,  the  climate  control  system,  etc.  Since  the  pressures  in 
aircraft  fuel  lines  generally  exceed  the  critical  pressure  of  the  fuel,  the  above 
mentioned  thermal  loads  will  cause  the  fuel  temperature  to  increase  beyond  the 
critical  temperature,  i.e.,  the  fuel  will  be  in  the  supercritical  regime. 

The  transport  and  chemical  behavior  of  fuels  under  supercritical  conditions  is  not 
well  understood.  Properties  such  as  density,  viscosity,  conductivity,  diffusivity, 
specific  heat,  etc.  undergo  steep  variations  near  the  critical  point.  These  variations 
are  expected  to  alter  the  heat  transfer  characteristics  of  the  fuel.  The  thermal  stability 
of  supercritical  fuels  is  of  great  concern.  If  the  temperature  of  the  fuel  exceeds  a 
certain  value,  the  fuel  undergoes  chemical  decomposition  to  form  insolubles  such 
as  gums  and  particle  clusters  that  get  deposited  on  fuel  lines,  valves,  etc.  Therefore, 
the  heat  management  has  to  be  done  in  a  temperature  range  over  which  the  fuel  is 
thermally  stable.  For  JP-8,  this  range  is  below  325°F.  Advanced  fuels  such  as  JP-900 
are  expected  to  be  thermally  stable  below  900®F. 
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Effective  heat  management  in  high  performance  aircraft  requires  a  better 
understanding  of  the  supercritical  transport  and  chemical  properties  of  the  fuel. 
The  Phase  I  study  addressed  the  issue  of  supercritical  heat  transfer.  An 
experimental /computational  approach  was  adopted  during  Phase  I.  The  Phase  n 
study  will  involve  the  analysis  of  precombustion  chemistry  and  deposit  formation 
under  supercritical  conditions. 

1.2  Phase  I  Objectives 

The  specific  objectives  of  the  Phase  I  study  were: 

a.  to  formulate  models  to  calculate  supercritical  transport  properties  such 
as  density,  viscosity,  conductivity  and  specific  heat; 

b.  to  conduct  experiments  with  Sulfur  Hexafluoride  (SF6)  as  a  surrogate 
fuel  and  to  measure  density  distributions  inside  the  test  cell; 

c.  to  obtain  functional  relationships  between  non-dimensional 
parameters  such  as  the  Nusselt  number  and  the  Re)molds  number; 

d.  to  perform  preliminary  validation  of  the  supercritical  transport 
models;  and 

e.  to  develop  recommendations  for  the  Phase  n  study. 

1.3  Phase  I  Accomplishments 

All  of  the  above  mentioned  objectives  were  accomplished.  The  Phase  I  work 
resulted  in  the  formulation  of  a  comprehensive  supercritical  transport  model.  This 
model  was  incorporated  into  an  existing  Computational  Fluid  Dynamics  (CFD)  code. 
This  has  resulted  in  a  unique  capability  to  compute  supercritical  flow  behavior  based 
on  local  conditions.  The  modified  code  was  used  to  simulate  heat  transfer  characteris¬ 
tics  over  a  range  of  conditions.  The  results  were  compared  with  data  obtained  from 
parallel  experiments  at  the  University  of  Iowa.  SFe  was  used  as  a  surrogate  fuel  for 
the  Phase  I  study. 

The  specific  achievements  of  the  Phase  I  study  are  listed  below. 
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a.  A  model  was  formulated  to  calculate  the  supercritical  transport 
properties  based  on  the  local  pressure,  temperature  and  composition. 
The  low  pressure  properties  were  calculated  from  the  kinetic  theory  of 
gases  and  appropriate  dense  gas  corrections  were  used  to  estimate  the 
properties  at  supercritical  conditions. 

b.  The  flow  inside  the  test  cell  was  experimentally  visualized  and 
shadowgraphs  were  obtained  to  show  density  variations. 

c.  Preliminary  validation  of  the  model  predictions  was  done  by 
comparing  the  results  against  existing  data  for  SF6.^ 

d.  The  supercritical  transport  model  was  incorporated  into  an  existing 
CFD  code,  REFLEQS  (REactive  FLow  EQuation  Solver).  The  code  was 
used  to  compute  heat  transfer  characteristics  in  a  2-D  channel.  The 
results  of  the  supercritical  model  were  compared  with  ideal  gas 
predictions  in  order  to  estimate  the  error  involved  with  ideal  gas 
approximations. 

e.  The  wall  heat  flux  was  computed  for  a  range  of  flow  conditions  in 
order  to  obtain  a  correlation  between  the  Nusselt  number  and  the 
Reynolds  number.  The  computed  relationship  was  compared  with 
standard  correlations  for  wall  heat  transfer. 

f.  A  series  of  parametric  simulations  were  performed  with  different 
gravity  vectors  to  analyze  the  effect  of  natural  convection.  This  study 
was  used  to  obtain  the  most  optimum  orientation  for  the  test  cell  in 
order  to  eliminate  flow  instabilities  generated  by  buoyancy. 

g.  A  3-D  model  was  used  to  simulate  the  flow  pattern  and  supercritical 
transport  inside  the  experimental  test  cell.  The  results  of  the 
simulations  were  compared  with  experimental  data. 

h.  A  preliminary  work  plan  for  Phase  n  was  prepared. 
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The  details  of  each  of  these  tasks  are  discussed  in  the  following  sections. 

1.4  Outline  of  the  Report 

The  next  section  (Section  2)  outlines  the  detailed  formulation  of  the  supercritical 
model.  Section  3  presents  the  computation  of  properties  for  SFg  at  various 
temperatures  using  the  model.  Computations  for  supercritical  heat  transfer  in  a  2-D 
channel  are  also  presented.  Section  4  presents  the  results  of  2-D  simulations  made 
to  estimate  the  effects  of  natural  convection  on  the  flow  patterns  inside  the  test  cell. 
Section  5  describes  the  experimental  setup  and  apparatus  and  the  data  obtained  from 
the  experiments.  Section  6  presents  the  results  of  3-D  simulations  of  the  test  cell  and 
comparisons  with  data.  The  conclusions  from  this  study  and  recommendations  for 
future  work  are  discussed  in  Section  7. 
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2. 


FORMULATION  OF  SUPERCRITICAL  TRANSPORT  MODEL 


This  section  presents  the  formulation  of  models  to  calculate  the  variation  of 
properties  (such  as  viscosity,  conductivity,  density,  specific  heat,  etc.)  in  the  critical 
and  supercritical  regimes. 

2.1  Governing  Equations 

The  mass  conservation  for  the  flow  can  be  represented  by  the  continuity  equation. 
This  formulation  accounts  for  a  variable  density  flowfield: 


^  +  V-pu  =  0.0 


(2.1) 


where  p  is  the  density  of  the  fluid  and  u  is  the  velocity  vector. 

The  momentum  conservation  equation  for  the  flow  is  written  as: 


+  Vpuu  =  -Vp  +  V-  Ti  (Vu  +  Vu^  +  Pg 


(2.2) 


where  t]  is  the  viscosity  and  g  is  the  gravity  vector.  It  is  assumed  that  the  gas 
behaves  as  a  Newtonian  fluid.  The  energy  conservation  equation  for  a  single 
component  system  (assuming  low  Mach  munber  and  neglecting  dissipation)  is 
given  by: 


(2.3) 


where  Cp,  A,  and  T  are  the  specific  heat,  thermal  conductivity,  and  temperature  of 
the  fluid,  respectively. 
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The  properties  such  as  p,  t],  X,  and  Cp  in  the  above  equations  have  to  be  calculated 
for  the  supercritical  regime.  Ideal  gas  approximations  are  not  valid  due  to  steep 
variations  (in  these  properties)  across  the  critical  point.  The  following  sections 
present  a  detailed  formulation  for  the  estimation  of  supercritical  properties. 

2.2  Calculation  of  Viscosity 

The  low  pressure  viscosity  is  calculated  from  the  expression  by  Chimg,  et  al^, 


T1 


3/2 


=  40.785 


vfa^ 


(2.4) 


7]  =  viscosity,  ^ 

M  =  Molecular  Weight,  g/mol 

T  =  Temperature,  °K 

Vc  =  critical  volume,  cm3 /mol 

Qg  is  di  collision  integral  of  the  form: 


a=  [a  {'P)'®]  +  C  [exp  (-DP)]  +  E[exp  {-FT)]  (2.5) 

where 


T  =  kT/e,  A  =  1.16145,  B  =  0.14874 
C  =  0.52487,  D  =  0.77320,  E  =  2.16178,  F  =  2.43787 
k  =  Boltzmann  Constant 

e  =  characteristic  energy  of  interaction  between  molecules 

The  collision  integral  is  based  on  a  Lennard-Jones  potential  function  of  the  form. 


\|/  (r)  =  4e 


(2.6) 
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where  a  is  the  collision  diameter  of  the  molecule. 


The  factor  elk  can  be  expressed  as: 


—  = - - —  (2  7) 

k  1.2593  ^ 

where  Tc  is  the  critical  temperature.  Therefore, 

T-*  =  1.2593  Tr  (2.8) 

where  Tr  =  T/Tc  is  the  reduced  temperature. 

The  factor  Fc  is  calculated  as 

f  ^  =  1  -  0.2756  0)  +  0.059035  \4  +  k  (2.9) 

where  cd  is  an  acentric  factor  dependent  on  the  polarity  of  the  molecule,  /Xr  is  a 
dimensionless  dipole  moment  and  fc  is  a  correction  for  highly  polar  molecules. 

The  high  pressure  viscosity  using  the  Chung  formulation  is  given  by: 


,363U[Mr) 

T1  =T1  - i ^ 


1/2 


u2/3 

c 


(2.10) 


where 


-.{ll 


112 


+  T\ 
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V  =  molar  volume  {cm^lmot\ 


Gi 


1  -05  y 

{1-y? 


G2  = 


E,  |[2  -  “P  +  £2  G^xp  (£5  y)  +  £3  G, 


Ej  E^  +  E2  +  E3 


T]**  =  Ejy^G2exp 


Eg  +  Eg 


(T] 


-1 


The  constants  Ei  through  Eio  are  tabulated  in  the  following  table. 


Table  2-1.  Constants  for  the  Viscosity  Model 


£,-  =  af  +  b^w  +  CfUi  +  dik 


/ 

a/ 

bi 

C/ 

di 

1 

6.234 

50.412 

-51.680 

1189.0 

2 

1.210  X  10*3 

-1 .154  X  10*3 

-6.257  X  10*3 

0.03728 

3 

5.283 

254.209 

-168.48 

3898.0 

4 

6.623 

38.096 

-8.464 

31.42 

5 

19.745 

7.630 

-14.354 

31.53 

6 

-1.900 

-12.537 

4.985 

-18.15 

7 

24.275 

3.450 

-11.291 

69.35 

8 

0.7972 

1.117 

0.01235 

-4.117 

9 

-0.2382 

0.06770 

-0.8163 

4.025 

10 

0.06863 

0.3479 

0.5926 

-0.727 

4240ri  t1 


The  above  formulation  for  calculating  dense  fluid  viscosities  shows  good  agreement 
with  experimental  data  and  the  errors  are  usually  less  than  5%. 
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The  high  pressure  thermal  conductivity  of  a  substance  is  given  by  the  expression  of 
Chung  et  al2 

A  =  thermal  conductivity  ,  W/(m-°K) 

Tjo  =  Low  pressure  gas  viscosity,  N  •  S/m2 
M'=  Molecular  weight,  kg/ mol 

q  =  3386  X  20'^  (^^')  /^2/3 
The  factor  f'  is  given  as: 

^=^1+ a  {j().215  +  028288a  - 1 .061  fi + 0.26665Z]I  [0.6366  +  )3Z+ 1 .061  a0\j  (2.12) 

where 

C.  3 
“  R'2 

p  =  0.7362  -  0.7109W  + 1.3168(0^ 
z  =  2.0  +  103T^^ 

Cv  =  Specific  heat  at  constant  volume 
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R  =  Universal  gas  constant 


(y )  (-B^y)]  +  ^2^1  exp  (B5y)  +  B3  Gj 

^2=  Bj  B4  +  B2  +  B3 

the  coefficients  Bi  to  By  are  tabulated  below. 

Table  2-2.  Constants  for  the  Thermal  Conductivity  Model 


/ 

a/ 

bi 

C/ 

d. 

1 

2.4166  E+0 

7.4824  E-1 

-9.1858  E-1 

1.2172  E+2 

2 

-5.0924  E-1 

-1.5094  E+0 

-4.9991  E+1 

6.9983  E+1 

3 

6.6107  e+0 

5.6207  E+0 

6.4760  E+1 

2.7039  E+1 

4 

1.4543  E+1 

-8.9139  E+0 

-5.6379  E+1 

7.4344  E+1 

5 

7.9274  E-1 

8.2019  E-1 

-6.9369  E-1 

6.3173  E+1 

6 

-5.8634  E+0 

1.2801  E+1 

9.5893  E+0 

6.5529  E+1 

7 

9.1089  E+1 

1.2811  E+2 

-5.4217  E+1 

5.2381  E+2 

4240/1 12 


2.4  Calculation  of  Density 


The  Corresponding  States  principle  is  used  to  calculate  the  density  of  the  fluid. 


Z 


-EY. 

~  RT 


where  Z  is  the  compressibility  factor.  Z  is  obtained  as: 

Z  =  (t,  ,Pr)  + (a  (t^  Pr) 


(2.13) 


(2.14) 
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The  function  Z(0)  would  apply  to  spherical  molecules  and  is  a  deviation 
function.  Z(0)  and  are  tabulated  as  a  function  of  T,  and  Pr  in  Tables  3-2  and  3-3  of 
Reference  2,  respectively. 

2.5  Calctdation  of  Enthalpy 

For  a  pure  substance. 


where 

h  =  enthalpy 
s  =  entropy 
V  -  specific  volume 

Also, 


dh  -  Tds  +  V  dp 


(2.15) 


(2.16) 


where 


Cp  =  Specific  heat  at  constant  pressure 


Substituting  Equation  (2.15)  in  Equation  (2.16)3, 


The  enthalpy  change  at  constant  temperature  is  given  as: 


(2.17) 


dh'p  — 


(2.18) 


ZRT 

Substituting  v  = - 


we  get 
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„  -RT^fdZ]  . 


=  -^T'Jr  ^ 


(2.19) 


The  zero  limit  for  pressure  signifies  the  enthalpy  of  an  ideal  gas,  h*.  The  real  gas 
enthalpy  ‘h’  at  any  temperature  can  be  obtained  from  the  above  relationship.  The 
RHS  of  Equation  (2.19)  is  a  function  of  T,.  and  Pf,  i-e. 


(2.20) 


where  F<0)  and  FC^)  are  tabulated  as  ftmctions  of  Tr  and  Pr  in  Tables  5-2  and  5-3  of 
Reference  2,  respectively.  The  ideal  gas  enthalpy  ft*  is  obtained  from  the  JANNAF 
data. 

It  is  observed  that  the  dependence  of  'h'  on  the  temperature  is  highly  non-linear. 
An  iterative  process  is  adopted  to  obtain  the  temperature  from  the  enthalpy  at  any 
given  point  in  the  computational  domain. 


Similar  to  enthalpy,  the  specific  heat  can  be  expressed  as: 


(2.21) 
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where  Cp*  is  the  ideal  gas  specific  heat.  Cp*  can  be  calculated  from  JANNAF  data. 
(ACp*)<0)  and  (ACp)(i>  are  tabulated  in  Tables  5-8  and  5-9  of  Reference  2,  respectively. 
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3. 


COMPUTATION  OF  SUPERCRITICAL  TRANSPORT  FOR  SF6 


Sulfur  Hexafluoride  (SF6)  was  selected  as  the  surrogate  fuel  for  the  Phase  I 
simulations  and  experiments.  SF6  is  non-toxic  and  non-flammable,  and  can  be 
easily  stored  as  a  liquid.  On  heating,  it  is  a  very  stable  compound  with  inert 
behavior  similar  to  nitrogen.  Also,  its  critical  state  can  be  achieved  easily  in  the 
laboratory.  This  makes  SF6  very  suitable  for  the  preliminary  experiments. 

The  critical  properties  for  SF6  are  : 

Tc  =  318.7°K 
Pc  =  37.6  bars 

Pc  =  734.7  kg/m3 
Zc  =  0.282 


The  acentric  factor  for  SF^  is  =  0.286  and  the  molecular  weight  of  SFg  is  146.05 
kg/mol. 

3.1  Property  Variation  with  Temperature 

The  property  calculations  were  done  at  the  critical  pressure  of  37.6  bars.  Figure  3-1 
shows  the  density  variation  with  temperature.  It  is  observed  that  the  density 
decreases  by  a  factor  of  3  across  the  critical  point.  The  symbols  represent  data  for  SFe^ 
in  the  subcritical  regime.  Reliable  data  are  not  available  in  the  critical  and 
supercritical  regimes.  The  comparison  shows  that  the  model  is  able  to  predict  the 
density  and  its  variation  reasonably  well  in  the  subcritical  regime.  Figure  3-2  shows 
the  dynamic  viscosity  variation  with  temperature.  A  steep  decrease  in  viscosity  is 
observed  near  the  critical  point.  This  correlates  well  with  other  observations.L4,5 
The  thermal  conductivity  and  specific  heat  variations  are  shown  in  Figures  3-3  and 
3-4,  respectively.  It  is  observed  that  both  quantities  have  peaks  in  the  vicinity  of  the 
critical  point.  Recent  investigations^  have  reported  this  unusual  behavior  for  the 
conductivity  and  specific  heat. 
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Figure  3-1.  Variation  of  Density  with  Temperature.  Symbols  represent  data  from 
Reference  1 


Figure  3-2.  Variation  of  Dynamic  Viscosity  with  Temperature 
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Specific  Heat  (J/kg  -  °K) 


Numerical  simulations  were  performed  to  calculate  supercritical  heat  transfer  in  a 
two-dimensional,  rectangular  channel.  The  schematic  of  the  channel  is  shown  in 
Figure  3-5a.  The  dimensions  of  the  channel  correspond  with  those  of  the  test  cell 
used  in  the  experiments.  The  flow  enters  at  a  subcritical  temperature  of  300®K  with 
a  velocity  of  1.16  mm/sec  (15cc/min).  The  bottom  wall  is  maintained  at  a 
supercritical  temperature  of  323°K.  Therefore,  the  wall  heating  causes  the  fluid  to 
become  critical  as  it  flows  along  the  wall.  The  computational  grid  (100  x  50)  is 
shown  in  Figure  3-5b. 


Figure  3-5a.  Schematic  of  the  Channel 


Figure  3-5b.  Computational  Grid  (100  x  50)  for  the  Channel 
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Figure  3-6  shows  the  temperature  distribution  in  the  channel.  Figure  3-7  shows  the 
density  distribution  in  the  channel.  It  is  observed  that  the  density  decreases  sharply 
near  the  wall  due  to  fluid  becoming  supercritical.  The  heat  flux  was  computed 
along  the  wall  and  is  shown  by  curve  1  in  Figure  3-8.  The  heat  flux  is  large  at  the 
leading  edge  of  the  wall  (due  to  thinner  botmdary  layers  at  the  leading  edge). 


temp  contours 

FMIN  3  a00E-i-02 

FMAX  3 . 2E8E-Ha2 

CONTOUR  LEVELS 

1  3  000E+02 

2  3  012E+02 

3  3  024E+02 

4  3 , 03SE-^02 

5  3  048E+02 

S  3 . 060E+02 

7  3.072E+02 

8  3 . 084E+02 

9  3 . 0S6E+02 

la  3 . 108E+02 

11  3.120E+02 

12  3 . 132E+02 

13  3 , 144E+02 

14  3  ,  15SE-H02 

15  3 . isaE+a2 

16  3 . 180E+02 

17  3  ,  192E+a2 

18  3 . 204E+«2 

19  3,2ieE-Ha2 

20  3  228E+a2 
QK> 


Figure  3-6.  Temperature  Distribution  (°K)  in  the  Channel 


RHO  CONTOURS 
FMIN  4,316E+«2 
FMAX  1 . 454E+03 
CONTOUR  LEVELS 

1  4.3 16E-MJ2 

2  4 . 854E+02 

3  5,392E+02 

4  5 . 93ae-M82 

5  6.469e+02 

S  7 . 007E+02 

7  7 .  S45E-M82 

8  8 . 083E+02 

9  9.621E+02 

10  9.159E+02 

11  9.S97E+02 

12  1 . a24E+03 

13  1 . 077E+<}3 

14  1.13 18+03 

15  1 . 185E+03 

16  1 . 23^+03 

17  1 . 293E+03 

18  1 . 34SE+03 

19  1 . 400E+03 

20  1 . 4548+03 

0K> 


Figure  3-7.  Density  (kg/m3)  Distribution  in  the  Channel 
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Figure  3-8.  Heat  Flux  Along  the  Wall  of  the  Channel:  Curve  1  -  Supercritical 
Model;  Curve  2  -  Ideal  Gas  Model 

3.3  Comparison  with  Ideal  GasiModel  Predictions 

An  ideal  gas  model  was  used  to  calculate  the  heat  transfer  in  the  2-D  chaimel  in 
order  to  estimate  the  error  involved  with  the  ideal  gas  approximation.  The  channel 
pressure  was  increased  to  248  bars  in  order  to  ensure  that  the  density  of  the  gas 
corresponded  to  the  density  of  the  subcritical  fluid  in  the  previous  case.  The 
velocity  is  maintained  at  1.16  mm /sec  so  that  the  mass  flow  rates  through  the 
channel  are  the  same  for  the  ideal  gas  model  and  the  supercritical  model.  The 
computed  wall  heat  flux  along  the  wall  from  the  ideal  gas  model  is  shown  by  curve 
2  in  Figure  3-8.  It  is  observed  that  the  supercritical  model  predicts  larger  heat  flux 
values  at  the  wall  for  the  same  boundary  conditions.  Figure  3-9  shows  the  density 
distribution  in  the  channel  predicted  by  the  ideal  gas  model.  It  is  observed  that  the 
ideal  gas  model  predicts  only  a  small  decrease  in  density  compared  to  the 
supercritical  model  (Figure  3-7). 
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RHO 

CONTOURS 

FMIN 

1 

352EH-03 

FMAX 

1 

454E+03 

CONTOUR 

LEVELS 

1 

1 

352E+03 

2 

1 

363E+03 

3 

1 

375E+03 

4 

1 

386E+03 

c 

1 

3  97  E +03 

6 

1 

409E+03 

7 

1 

'‘20E+03 

a 

1 

43 lE+03 

9 

1 

443E+03 

10 

1 

454E+03 

0K> 


Figure  3-9.  Density  (kg/m3)  Distribution  (with  Ideal  Gas  Model)  in  the  Channel 
3.4  Correlation  Between  Dimensionless  Quantities 

A  series  of  simulations  were  performed  for  different  flow  rates  through  a  2-D 
channel  in  order  to  correlate  the  Nusselt  number  with  the  Reynolds  number.  The 
schematic  of  the  channel  is  shown  in  Figure  3-10.  A  fully  developed  velocity  profile 
is  assumed  at  the  inlet.  The  fluid  enters  at  a  temperature  of  300®K  and  the  walls  of 
the  channel  are  maintained  at  323®K.  The  mean  velocity  through  the  channel  is 
used  as  a  parameter  to  vary  the  Reynolds  number.  The  computed  wall  heat  flux  (for 
the  entire  wall)  is  plotted  (curve  1)  in  Figure  3-11  as  a  function  of  the  mean  velocity. 
Curve  2  shows  the  results  obtained  from  the  ideal  gas  approximation. 


Mean  Velocity 
=  mnVsec 

T  =  300‘’K 


Figure  3-10.  Schematic  of  the  Charmel  and  Boundary  Conditions 
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1  0e-i-04^ 


Vm  ,  Mean  Velocity  of  Flow  (mm/sec) 


Figure  3-11.  Variation  of  Heat  Flux  with  Mean  Velocity  Through  the  Channel. 
Curve  1  -  Supercritical  Model;  Curve  2  -  Ideal  Gas  Model 

A  curve  fit  of  the  heat  flux  values  results  in  the  following  correlations: 

=  (3.1) 

for  the  supercritical  model,  and 

(3.2) 

for  the  ideal  gas  model.  Where  is  the  heat  flux  over  the  length  L,  Vm  is  the  mean 
velocity,  and  Kj  and  Kz  are  proportionality  constants.  It  is  interesting  to  note  that 
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the  supercritical  model  and  the  ideal  gas  model  predict  the  same  dependence  for 
heat  flux  on  the  mean  velocity.  However,  the  proportionality  constant,  Kj,  for  the 
supercritical  model  is  larger  than  Kz,  the  proportionality  constant  for  the  ideal  gas 
model. 

The  Nusselt  number  is  expressed  as: 

=  T  (3.3) 

U 

where  h  is  the  heat  transfer  coefficient  =  ^  ,  L  is  the  length  of  the  wall,  and  X  is 

the  thermal  conductivity  of  the  freestream  fluid.  The  Reynolds  number  is  expressed 
as: 


(3.4) 


where  H  is  the  channel  height,  and  p  and  "n  are  the  density  and  dynamic  viscosity  of 
the  freestream  fluid.  Therefore,  Equations  (3.1)  and  (3.2)  can  be  rewritten  as: 


(3.5) 


The  standard  correlation  for  heat  transfer  in  a  long  duct  is  expressed  as:^ 


Nui  =  1.86lRe- Pr  - —\ 

where  Pr  is  the  Prandtl  number,  fib  and  fis  sie  the  dynamic  viscosities  in  the  bulk 
fluid  and  heated  surface,  respectively.  It  is  observed  that  the  supercritical  model  and 


pj  (3.6) 
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the  ideal  gas  model  predict  correctly  the  dependence  of  Nusselt  number  on  the 
Reynolds  number.  The  proportionality  constant  for  the  supercritical  model  is 

higher  than  that  of  the  ideal  gas  model.  From  Figure  3-11,  it  is  observed  that  Ki  == 
6K2.  However,  the  thermal  conductivity  predicted  by  the  supercritical  model  at  the 
heated  wall  is  only  3  times  as  large  as  that  predicted  by  the  ideal  gas  model.  In  order 
to  understand  the  sixfold  increase  in  heat  flux,  it  is  necessary  to  study  the  boimdary 
layer  profiles  in  the  near  wall  region.  Figure  3-12  shows  the  velocity  profiles 
(midway  in  the  channel)  predicted  by  the  two  models.  It  is  observed  that  the  steep 
decrease  in  density  predicted  by  the  supercritical  model  causes  the  fluid  near  the 
wall  to  accelerate  leading  to  a  thinner  boundary  layer  (reflected  by  the  larger  slope  in 
the  profile  near  the  wall).  This  results  in  a  larger  heat  flux  than  what  the  ratio  of 
conductivities  would  indicate. 


Velocity  (mm/sec) 


Figure  3-12.  Velocity  Profile  in  the  Channel.  - Ideal  Gas  Model;  —  Supercritical 

Model 
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4. 


ASSESSMENT  OF  NATURAL  CONVECTION  EFFECTS 


The  results  obtained  in  the  previous  section  show  that  the  density  of  SF6  decreases 
by  a  factor  of  3  across  the  critical  point.  This  decrease  occurs  over  a  temperature 
interval  of  10°K.  One  of  the  main  concerns  in  experimenting  with  supercritical 
fluids  has  been  the  effect  of  buoyancy  which  may  be  quite  significant  due  to  the 
presence  of  the  large  density  gradient.  Buoyancy  could  introduce  gravity  driven 
flow  instabilities  and  thus  complicate  the  study  (and  hence  the  understanding)  of 
supercritical  transport  phenomena.  One  of  the  goals  of  this  study  was  to  obtain  an 
appropriate  orientation  for  the  test  cell  that  would  minimize  flow  instabilities 
(especially  in  the  near  wall  region)  due  to  natural  convection.  This  section  presents 
the  results  for  flow  in  a  2-D  channel  (dimensions  and  boundary  conditions  are 
shown  in  Figure  3-5)  with  different  orientations  of  the  gravity  vector. 

4.1  Case  1 ;  Gravity  Vector  gy  =  -9.81  m/sec2 


Figure  4-1  shows  the  pressure  distribution  in  the  chaimel.  The  pressure  increases  in 
the  negative  y-direction  due  to  the  hydrostatic  pressure  gradient.  Figure  4-2  shows 
the  streamlines  for  this  case.  It  is  observed  that  the  pressure  gradient  interacts  with 
the  density  gradient  near  the  wall  to  produce  gravity  driven  recirculation  zones  (i.e., 
flow  separation  from  the  wall)  in  the  channel.  Figure  4-3  shows  the  heat  flux 
distribution  along  the  channel  wall.  The  recirculating  flow  results  in  a  highly  non¬ 
smooth  distribution  of  heat  flux  along  the  channel.  The  result  from  the  zero  gravity 
case  is  shown  as  a  dashed  line. 


PRES  CONTOURS 
FMIN  -3.724E+«0 

FMAX  I . 355E+«a 

CONTOUR  LEVELS 
1  -9.724e-»«0 

Z  1.308E-)^1 

3  3.5aaE+«i 

4  S.868E-^«1 

5  a.i4aE+«i 

6  1 . 043E+«2 

7  I.Z7  lE-naa 

8  1 . 49SE+02 

9  1 . 727E+02 

10  1 . 955E-Ma2 
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Figure  4-1.  Pressure  Distribution  (Pascals)  in  the  Channel  for  Case  1 
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9 

-3  363E-01 

10 

-2.820E-01 

11 

-2.277E-01 

12 

-1.73SE-01 

13 

-1.  192E-01 

14 

-S  498E-0a 

15 

-1.072E-02 

IE 

4 . 354E-02 

17 

3 . 780e-02 

18 

1.521E-01 

19 

2.063E-01 

20 

0K> 

2.606E-01 

Figure  4-2.  Streamline  Contours  for  Case  1 


1,  Length  Along  Channel,  (m) 

Figure  4-3.  Heat  Flux  Along  the  Channel  Wall.  -  Results  for  Case  1 

Results  with  Zero  Gravity 


4.2  Case  2  ;  Gravity  Vector  gy  =  9.81  m/serf 


The  direction  of  the  hydrostatic  pressure  gradient  is  reversed  as  compared  to  the 
previous  case.  Figures  4-4  and  4-5  show  the  pressure  distribution  and  the 
streamlines  in  the  channel  respectively.  Again,  flow  instabilities  due  to  buoyancy 
are  observed  in  the  channel.  The  effect  of  the  instabilities  are  reflected  in  the  heat 
flux  curve  shown  in  Figure  4-6. 


PRES  CONTOURS 
PMIN  -2.036E+«2 
FMAX  9.033E-01 
CONTOUR  LEVELS 

1  -2 . 036E+02 

2  - 1 . 8C9E+02 

3  - 1 . 5a2£-H32 

4  -1.3S4E-Ha2 

5  -1.127E+«2 
S  -a.998E+«l 
7  -S,72SE-Mai 
a  -4.454E+01 
S  -2.182E+01 

10  9.033E-0X 
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Figure  4-4.  Pressure  Distribution  (Pascals)  for  Case  2 
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Figure  4-5.  Streamline  Contours  for  Case  2 
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Figure  4-6.  Heat  Flux  Al  ..ng  the  Channel  Wall.  -  Results  for  Case  2; 

—  Results  with  Zero  Gravity 

4.3  Case  3 ;  Gravity  Vector  =  9.81  m/sec2 


The  gravity  vector  now  points  in  the  positive  x-direction,  i.e.,  the  channel  is  vertical 


with  the  flow  going  from  top  to  bottom.  Figure  4-7  shows  the  pressure  distribution 
in  the  channel.  The  streamlines  shown  in  Figure  4-8  indicate  the  presence  of 


buoyancy  driven  cells.  Figure  4-9  shows  the  heat  flux  distribution 


channel. 
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Figure  4-7.  Pressure  Distribution  (Pascals)  for  Case  3 
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Figure  4-8.  Streamline  Contours  for  Case  3 


1,  Length  Along  Channel,  (m) 


Figure  4-9.  Heat  Flux  Along  the  Channel  Wall. 
Results  with  Zero  Gravity 


Results  for  Case  3, 


4.4  Case  4  :  Gravity  Vector  gg  =  -9.81  m/sec2 


The  gravity  vector  is  in  the  negative  x-direction,  i.e.,  the  channel  is  vertical  with  the 
flow  going  from  the  bottom  to  the  top.  Figure  4-10  shows  the  pressure  distribution 
in  the  channel.  The  streamlines  in  Figure  4-11  indicate  a  smooth  flow  pattern  (i.e., 
there  is  no  flow  separation).  The  effect  of  buoyancy  in  this  orientation  is  to 
accelerate  the  low  density  fluid  near  the  wall  upwards.  However,  the  recirculation 
zones  observed  in  all  the  previous  cases,  are  absent.  This  is  also  reflected  in  the 
smooth  nature  of  the  curve  for  the  heat  flux  distribution  (Figure  4-12).  The  effect  of 
flow  acceleration  near  the  wall  produces  higher  values  for  the  heat  flux  compared  to 
the  zero  gravity  case. 
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Figure  4-10.  Pressure  Distribution  (Pascals)  for  Case  4 


STW1  CONTOliS 
PTIIN  S.199e-94 

FMAX  1.794e-91 

CQHTOfl  LEVELS 

1  S.  1S9E-94 

2  g.S53E-«3 

3  1.849e-92 

4  2.742E-92 

5  3.63^-92 

S  4.52Se-92 

7  5.422E-92 

a  S.313E-92 

9  7 . 299e-92 

19  8 .  192E-92 

11  8.39GE-92 

12  9.889E-92 

13  1.978e-«l 

14  l.lBaE-9l 

15  1.2S7e-91 

16  1.346E-91 

17  1.436E-91 

18  1.S2K-91 

13  1.S14E-91 

29  1.794E-91 

0K> 


Figure  4-11.  Streamline  Contoirrs  for  Case  4 
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1,  Length  Along  Channel  (m) 

Figure  4-12.  Heat  Flux  Along  the  Channel  Wall.  -  Results  for  Case  4;  — 

Results  with  Zero  Gravity 

From  the  above  simulations,  it  is  clear  that  the  best  possible  orientation  is  that 
corresponding  to  Case  4,  i.e.,  a  vertical  channel  with  flow  going  from  bottom  to  top. 
The  mean  flow  patterns  for  this  case  are  smooth  (without  flow  separation)  which 
result  in  a  smooth  heat  flux  distribution.  Again,  it  should  be  noted  that  the  effect  of 
gravity  in  this  orientation  is  to  accelerate  the  low  density  fluid  upwards  causing 
higher  rates  of  heat  transfer  from  (almost  an  order  of  magnitude  higher)  the  wall. 
Figure  4-13  shows  the  temperature  distribution  in  the  channel  for  case  4.  It  is 
observed  that  the  accelerating  flow  near  the  wall  causes  the  thermal  boundary  layer 
to  become  very  thin. 

Based  on  the  above  simulations,  a  vertical  orientation  (with  the  flow  going  from  the 
bottom  to  the  top)  was  selected.  The  experiments  discussed  in  the  next  section  were 
done  with  this  orientation  for  the  test  cell. 
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Figure  4-13.  Temperature  Distribution  (®K)  in  the  Channel  for  Case  4 
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5.  EXPERIMENTAL  RESULTS  FOR  SUPERCRITICAL  HEAT  TRANSFER 

This  section  presents  the  experimental  results  obtained  at  the  University  of  Iowa. 
SF^  was  used  as  the  simulant  fluid  and  the  experiments  were  performed  in  the 
subcritical,  critical,  and  supercritical  regimes.  The  results  are  shown  in  the  form  of 
shadowgraphs. 

5.1  Experimental  Apparatus 

A  square  charmel  with  an  interior  cross  section  of  15  mm  by  15  mm  is  used  for  the 
heat  transfer  experiments.  This  square  channel  is  fitted  with  removable  walls.  Two 
optical  quality  windows  are  fitted  on  the  opposite  sides  of  the  apparatus.  One  side  of 
the  channel  wall  is  fitted  with  a  heater  block,  and  the  other  side  a  thermocouple 
holder  block.  The  viewing  area  of  the  optical  access  is  15  mm  by  141  mm.  This 
I  square  channel  is  the  test  cell  for  the  heat  transfer  experiments.  The  heat  transfer 

test  cell  is  mounted  vertically  upward  on  a  vibration-isolated  optical  table.  A 
schematic  diagram  of  the  experimental  set-up  is  shown  in  Figure  5-1.  A  digitized 
(video)  image  of  the  apparatus  is  shown  in  Figure  5-2.  The  inlet  of  the  test  cell  is  its 
I  lower  port  and  the  exit  is  the  upper  one.  The  inlet  port  is  plumbed  to  a  feeder 

I  cylinder  which  is  essentially  a  high  pressure  canister.  The  interior  volume  of  this 

canister  is  2.5  liters.  The  feeder  cylinder  is  connected  to  liquid  SF6  supply  and  to  a 
high  pressure  helium  cylinder.  Liquid  SFe  is  supplied  to  the  feeder  cylinder  by  its 
own  vapor  pressure.  Helium  is  used  to  pressurize  the  liquid  SFg  in  the  feeder 
!  cylinder  to  desired  pressures.  A  pressure  relief  valve  set  at  4.34  MPa  (600  psia)  is 

installed  in-line  between  the  feeder  cylinder  and  the  helium  supply,  and  the  other 
installed  in  between  the  feeder  cylinder  and  the  test  cell.  The  pressure  is  measured 
by  a  pressure  transducer  (Omega  Model  PX302-1KGV)  and  monitored  by  a 
transducer  indicator  (Omega  Model  DP-350). 
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The  heater  block  is  electrically  heated  by  using  a  Variac  (Staco  Energy  Product,  Model 
3PN1010V).  K-type  thermocouples  are  spot-welded  onto  the  heater  block  to  monitor 
the  heater  block  temperature.  The  bulk  fluid  temperature  inside  the  test  cell  is 
measured  by  three  sheathed  K-type  thermocouples;  they  are  respectively  located 
near  the  inlet,  nud-section,  and  exit  of  the  test  cell.  These  sheathed  thermocouples 
are  introduced  to  the  flow  through  the  holder  block  which  is  located  at  the  opposite 
side  from  the  heater  block.  A  thermocouple  indicator  (Fluke  Model  21 90 A)  is  used 
to  monitor  the  temperatures  of  the  heater  block  and  the  bulk  fluid. 

The  SF6  flow  rate  is  regulated  by  a  fine  metering  valve  (Nupro  B-4MG),  and 
monitored  by  a  mass  flowrate  meter  (with  flow  elements  of  Hastings  Model  FST  and 
flowrate  monitor  Hastings  Model  PR-4AJ).  The  flow  meter  was  calibrated  for  air 
flow  rates.  The  conversion  factor  was  provided  by  the  vendor  for  SF6. 

The  shadowgraph  flow  visualization  is  used  for  the  present  study  of  the  transport 
process.  An  argon-ion  laser  (Coherent  Model  Innova-70)  is  used  as  the  light  source. 
The  optical  set-up  consists  of  a  condensing  lens  and  a  piano-cylindrical  lens.  The 
cylindrical  lens  is  used  to  accommodate  the  slender  rectangular  shape  of  the 
viewing  window.  The  shadowgraph  is  projected  onto  a  screen  and  is  recorded  by  a 
CCD  video  camera  (Panasonic  Model  AG-180).  A  frame  grabber  (Scion  Image 
Capture  H)  is  used  to  digitize  the  video  image  for  review  and  analysis.  A  still 
camera  is  also  used  to  record  some  of  the  shadowgraph  images. 

5.2  Results  and  Discussion 

At  zero,  or  low  SFg  flow  rates,  the  heating  of  SFe  at  a  near  critical  pressure  condition 
with  the  heater  block  maintained  at  the  near  critical  temperature  results  in 
recirculation  cells  in  the  channel.  At  a  zero  flow  rate,  a  stratified  region  is  first 
formed  in  the  near  exit  region  followed  by  a  gradual  down-flow  along  the  cold 
surface  of  the  unheated  sides.  This  down-flow  is  a  time-dependent  phenomenon, 
and  is  evident  in  the  video  taped  images  at  16.4  cc/min  inlet  flow  rate,  cf..  Figure  5- 
3a.  The  SF6  volumetric  flowrate  is  calculated  from  the  mass  flow  rate  with  the 
correction  for  the  compressibility  factor  difference  between  the  atmospheric 
condition  and  the  inlet  fluid  state. 
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Figure  5-3.  Shadowgraph  Images  of  the  Flow  in  the  Test  Cell 
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Vigorous  mixing  between  the  fluid  adjacent  to  the  heated  surface  and  the  "free- 
stream  fluid"  is  evident  from  the  ejection  phenomenon  of  the  fluid  parcels  (or 
blobs)  from  the  boundary  in  the  shadowgraph.  This  ejection  of  the  fluid  blobs  is 
similar  to  the  turbulence  spots  in  boundary-layer  flows.  Figure  5-3a  also  shows  the 
down  flow  along  the  unheated  surface.  The  down  flow  action  is  a  result  of  the 
cooling  of  the  hot  SF6  or  the  stratified  fluid  that  has  accumulated  in  the  upper 
region  of  the  channel.  Without  introducing  an  SFg  flow,  the  channel  flow  becomes 
a  cavity  flow  with  a  strong  mixing  of  the  fluid  inside  the  channel.  When  the  flow 
rate  is  introduced,  for  example,  at  16.4  cc/nun  of  Figure  5-3a,  the  lower  1/3  region  of 
the  field  of  view  shows  a  "calm  flow"  region.  This  "calm  flow"  region  grows  in  size 
when  the  SF6  flow  rate  is  increased  from  zero  imtil  a  linuting  condition  is  reached 
beyond  which  its  size  reduces  when  the  flow  rate  is  increased.  This  is  shown  in 
Figure  5-3b  where  the  flow  rate  at  the  inlet  is  maintained  at  35.6  cc/min.  At  this 
condition,  the  onset  of  the  "turbulence  spot"  in  the  boundary  layer  moves  toward 
the  channel  inlet.  With  further  increase  of  the  SFs  flow  rate,  the  onset  of  the  wavy 
boundary  layer  or  the  turbulence  spots  moves  toward  the  channel  inlet,  for  example 
as  shown  in  Figures  5-3c  and  5-3d  where  the  inlet  flow  rate  is  maintained  at  41.0  and 
38.3  cc/min.  The  system  pressure  is  maintained  at  Pr  =  1.02,  1.01,  ^d  0.83  for 
Figures  5-3b,  5-3c,  and  5-3d,  respectively;  where  Pr  is  the  reduced  pressure  based  on 
the  critical  pressure  of  SFe-  Figure  5-3e  shows  the  turbulence  spots  develop  into  a 
turbulent  boimdary-layer-like  flow  shown  by  the  large  structures  in  shadowgraph 
visucdization.  It  should  be  noted  that  the  ejection  of  fluid  parcels  is  also  observed 
along  the  unheated  surface  due  to  the  cooling  effect.  This  ejection  of  the  fluid  blobs, 
however,  may  be  an  illusion  of  the  collapsing  fluid  parcels  due  to  decrease  of  the 
specific  volume  when  the  hot  fluid  approaches  the  cold  surface.  When  the  heater 
block  is  maintained  at  40®C  (Tr  =  0.98),  and  the  flow  rate  maintained  at  30.1  cc/min, 
the  heating  increases  the  liquid  SF6  from  21®C  at  the  inlet  to  30®C  at  the  top  of  the 
heater  block. 

When  the  heater  block  is  maintained  at  a  high  temperature,  e.g.,  Figtu-e  5-3f,  where 
the  mid-section  thermocouple  registers  a  temperature  of  Tr  =  1.02  and  the  exit  fluid 
temperature  of  Tr  =  1,04,  the  opalescence,  a  near  critical  phenomenon,  is  seen.  Tr  is 
the  reduced  temperature  based  on  the  critical  temperature  of  SFe.  The  opalescence 
results  in  a  highly  distorted  image,  and  a  nearly  opaque  section  of  Figure  5-3f  due  to 
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the  large  variation  of  the  refractive  index  due  to  the  fluid  approaching  the  critical 
point.  The  flowrate  of  Figure  5-3f  is  maintained  at  30.1  cc/min.  This  large  variation 
of  the  refractive  index  is  evident  in  the  shadowgraph  of  the  near  wall  fluid 
structures  that  shows  the  refraction  of  the  light  beam  outside  the  boundary  of  the 
viewing  port. 
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6. 


3-D  COMPUTATIONS  OF  SUPERCRITICAL  TRANSPORT  IN  THE  TEST 
CELL 


The  supercritical  models,  described  in  Section  2,  were  incorporated  into  a  3-D  CFD 
code.  In  the  experiments,  the  fluid  enters  the  test  cell  (square  in  cross  section) 
through  a  cylindrical  tube.  This  gives  rise  to  a  3-dimensional  flow  near  the  inlet 
thus  necessitating  the  use  of  a  3-D  model. 

6.1  Computational  Grid  and  Boundary  Conditions 

The  schematic  of  the  test  cell  is  shown  in  Figure  6-1.  The  flow  enters  through  a 
circular  tube  into  the  rectangular  test  cell.  A  fully  developed  parabolic  velocity 
profile  is  specified  at  the  inlet  to  the  test  cell.  A  mean  velocity  of  13  mm /sec 
(corresponding  to  a  flow  rate  of  25  1/m)  is  used  as  the  inlet  flow  condition.  The  inlet 
temperature  is  300®K  and  the  bottom  wall  of  the  test  cell  is  maintained  at  323°K. 
The  test  cell  pressure  is  37.6  bars.  The  flow,  after  passing  through  the  test  cell,  exits 
into  a  circular  tube.  The  orientation  of  the  test  cell  (based  on  earlier  2-D 
simulations)  is  vertical  with  the  flow  going  from  the  bottom  to  the  top,  i.e.,  the 
gravity  vector  is  opposite  to  the  mean  flow  direction. 


Figure  6-1.  Schematic  of  the  Test  Cell 
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The  computational  grid  for  the  test  cell  is  shown  in  Figvne  6-2.  The  grid  has  61 
nodes  along  the  tube,  43  cells  in  the  direction  normal  to  the  heated  plate,  and  44  cells 
in  the  transverse  plane  of  the  heated  plate.  The  grid  is  clustered  near  the  heated 
boundary  in  order  to  resolve  steep  variations  in  the  profiles  of  the  flow  properties. 
A  time  dependent,  3-D  calculation  was  performed  to  analyze  the  transport  process 
within  the  test  cell. 

6.2  Results  of  Numerical  Simulations 

Figure  6-3  shows  the  density  distribution  in  the  test  cell.  It  is  observed  that  there  is  a 
steep  variation  in  density  near  the  heated  wall.  However,  this  variation  is  confined 
to  a  very  thin  layer  adjacent  to  the  wall.  The  velocity  vectors,  shown  in  Figure  6-4, 
indicate  the  flow  acceleration  near  the  heated  wall  due  to  buoyancy  effects.  The 
accelerating  flow  causes  the  boundary  layer  to  become  extremely  thin.  The 
impingement  of  this  high  velocity  flow  on  the  rear  wall  of  the  test  cell  results  in  a 
recirculation  zone  near  the  exit.  The  computed  flow  results  show  a  downflow  of 
fluid  on  the  unheated  wall.  This  is  also  observed  in  the  experiments  (Figure  5-3). 
Figure  6-5  shows  the  temperature  distribution  inside  the  test  cell.  The  increase  in 
temperature  is  confined  to  a  thin  layer  adjacent  to  the  wall.  However,  the 
recirculation  zone  near  the  exit  brings  the  heated  fluid  from  the  wall  to  the  core 
region  of  the  test  cell. 

Figure  6-6  shows  the  heat  flux  distribution  on  the  heated  wall  of  the  test  cell.  The 
heat  flux  values  are  higher  at  the  leading  edge  of  the  wall  due  to  thinner  boundary 
layers.  However,  the  general  profile  of  the  heat  flux  distribution  is  the  same  as  that 
predicted  by  the  2-D  simulations.  Some  3-dimensionalities  are  observed  in  the  heat 
flux  distribution  due  to  the  effect  of  the  velocity  profiles. 

6.3  Discussion 

The  major  features  of  the  experimented  flow  field  are  reproduced  by  the  model.  The 
acceleration  of  the  fluid  along  the  heated  wail,  the  downflow  along  the  unheated 
wall,  and  the  recirculating  flow  near  the  test  cell  exit  are  predicted  correctly  by  the 
simulations. 
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Figure  6-2.  Computational  Grid  for  the  Test  Cell 


6-3.  Density  Distribution  (kg/cu.m)  in  the  Test  Cell 
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Figure  6-4.  Velocity  Field  (m/sec)  inside  the  Test  Cell 


igure  6-5.  Temperature  Distribution  (degrees 
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It  is  observed  that  supercritical  phenomena  manifest  themselves  in  a  complex 
manner  in  the  flow  field.  It  is  difficult  to  understand  the  interactions  between 
various  physical  mechanisms  through  experiments  alone.  The  modeling  of  these 
flows  yields  considerable  insight  into  the  governing  physical  processes.  Experiments 
are  needed  to  judge  the  validity  of  the  model  and  to  guide  model  development. 
The  current  Phase  I  effort  has  demonstrated  proof-of-concept  that  joint 
computational /experimental  studies  are  crucial  in  order  to  develop  and  improve 
predictive  capabilities  for  supercritical  flows. 
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7. 


CONCLUSIONS  AND  RECOMMENDATIONS 


This  section  summarizes  the  work  performed,  outlines  the  main  findings  of  this 
study  and  presents  recommendations  for  future  work. 

7.1  Summary  of  Work  Performed 

All  of  the  objectives  of  the  Phase  I  study  have  been  met.  The  main 
accomplishments  of  the  project  are  summarized  below. 

a.  Development  of  Supercritical  Transport  Model:  A  comprehensive 
model  was  formulated  to  calculate  the  transport  properties  of  SF6  as  a 
function  of  local  pressure,  temperature  and  composition.  This  model 
was  able  to  reproduce  correctly  the  variation  of  these  properties  in  the 
vicinity  of  the  critical  point. 

b.  Incorporation  of  Model  into  2-D  CFD  Code:  The  supercritical  transport 
model  was  incorporated  into  a  general  purpose  2-D  CFD  code  which 
was  then  used  to  perform  heat  transfer  calculations  in  a  channel.  Wall 
heat  fluxes  were  computed  along  the  channel  wall  and  compared  with 
predictions  from  the  ideal  gas  model.  It  was  demonstrated  in  this  study 
that  the  ideal  gas  model  is  highly  inappropriate  for  modeling 
supercritical  flows  due  to  its  inability  to  resolve  the  steep  variations  in 
properties  around  the  critical  point.  The  model  was  also  used  to  study 
the  relationship  between  dimensionless  quantities  such  as  the  Nusselt 
number  and  the  Reynolds  number. 

c.  Assessment  of  Natural  Convection  Effects:  The  steep  variation  in 
density  aroimd  the  critical  point  gives  rise  to  strong  buoyancy  effects  in 
low  speed  supercritical  flows  imder  laboratory  conditions.  One  of  the 
objectives  of  this  study  was  to  run  the  experiments  at  a  low  enough 
Reynolds  number  (to  avoid  turbulent  flow  situations)  and  yet  prevent 
gravity  driven  flow  instabilities  from  occurring  within  the  test  cell. 
The  supercritical  model  was  used  to  obtain  the  test  cell  orientation  that 
best  met  these  objectives.  A  vertical  orientation  of  the  test  cell  with  the 
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flow  going  from  the  bottom  to  the  top  was  found  to  minimize 
buoyancy  driven  instabilities. 


d.  Experiments  with  SFg:  Experiments  were  carried  out  (as  part  of  a 
subcontract)  by  Professor  L.  D.  Chen  at  the  University  of  Iowa.  The 
flow  inside  the  test  cell  was  experimentally  visualized  in  the 
subcritical,  critical,  and  supercritical  regimes.  The  results  were 
presented  in  the  form  of  shadowgraphs.  The  experimental  data  were 
used  for  preliminary  validation  of  the  model. 

e.  Incorporation  of  Model  into  3*D  CFD  Code:  The  supercritical  fluid 
model  was  incorporated  into  a  3-D  CFD  code  to  simulate  the  flow  in 
the  test  cell.  Since  the  geometry  of  the  test  cell  produces  a  3- 
dimensional  flow,  a  3-D  model  was  necessary  for  the  simulations.  The 
results  of  the  simulations  were  compared  with  the  data  obtained  from 
experiments. 

The  Phase  I  effort  has  resulted  in  the  development  of  a  novel  and  unique  capability 
to  simulate  critical  and  supercritical  flows  using  a  general  purpose  CFD  code. 
Therefore,  the  effect  of  complex  boundary  and  operating  conditions  can  be 
simulated  without  the  use  of  adhoc  empirical  correlations. 

7.2  Main  Findings  from  the  Phase  I  Study 

Since  the  properties  of  supercritical  fluids  are  strong  functions  of  the  local  pressure, 
temperature  and  composition,  it  is  necessary  to  incorporate  supercritical  models 
into  CFD  codes  capable  of  resolving  the  interactions  between  these  quantities. 
Phenomena  such  as  heat  transfer  zmd  chenucal  reactions  are  an  outcome  of  these 
interactions.  The  results  of  this  study  show  that  heat  transfer  characteristics  are 
altered  by  an  order  of  magnitude  due  to  supercritical  effects.  This  behavior  is  certain 
to  affect  the  temperature  distribution  in  the  test  cell  and  hence  the  chemistry 
(through  the  non-linear  dependence  on  temperature).  It  was  also  found  that 
phenomena  such  as  gravity  should  be  given  careful  consideration  in  conjunction 
with  supercritical  fluids.  The  steep  density  gradients  occurring  in  supercritical  flows 
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may  give  rise  to  intrinsic  flow  instabilities  under  the  influence  of  pressure  gradients 
(hydrostatic  and  hydrodynamic). 

The  Phase  I  work  has  shown  the  proof-of-concept  of  coupling  detailed  supercritical 
fluid  models  with  conventional  CFD  codes  to  predict  flow  and  heat  transfer 
characteristics.  This  approach  provides  the  capability  of  analyzing  the  effects  of 
complex  boundary  conditions  on  local  flow  behavior.  It  also  provides  a  base  for 
including  more  physical  models  for  phenomena  such  as  chemical  reactions  in  order 
to  study  precombustion  chemistry  in  jet  fuels. 

7.3  Recommendations  for  Phase  II 

In  order  to  develop  a  general  purpose  capability  to  simulate  the  precombustion 
chemistry  and  transport  phenomena  in  supercritical  jet  fuels,  the  following 
refinements  are  being  considered  for  Phase  11  of  this  project: 

a.  Development  of  Multi-Component  Capability:  The  Phase  I  effort 
focused  on  the  analysis  of  the  supercritical  behavior  of  SF6  which  was 
used  as  a  surrogate  fuel  during  this  study.  Jet  fuels  on  the  other  hand 
are  a  complex  mixture  of  various  hydrocarbons.  This  requires  the 
capability  to  model  additional  features  such  as  multi-component 
diffusion,  thermodiffusion  (the  Soret  Effect),  the  Dufour  Effect,  etc.  in 
the  supercritical  regime. 

b.  Development  of  Chemistry  and  Deposition  Models:  The 

precombustion  chemistry  in  jet  fuels  involves  a  range  of  complex 
processes  such  as  (i)  auto-oxidation  of  the  fuel  at  low  temperatures 
with  the  dissolved  oxygen,  (ii)  thermal  pyrolysis  of  the  hydrocarbon 
molecules  at  higher  temperatures,  (iii)  formation  of  insolubles  such  as 
gums,  particles,  etc.,  (iv)  transport  and  deposition  of  the  insolubles  on 
the  walls  of  the  fuel  lines.  It  is  necessary  to  develop  models  for  the 
homogeneous  gas  phase  reactions  that  give  rise  to  the  insolubles  and 
the  heterogeneous  reactions  that  actually  result  in  surface  deposition. 
The  study  of  particle  formation  should  also  include  the  analysis  of 
homogeneous  and  heterogeneous  nucleation. 
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c.  Development  of  Turbxxlence  Models  for  Supercritical  Flows:  All  of  the 
modeling  work  in  Phase  I  was  done  using  laminar  flow  models.  From 
practical  considerations,  it  is  necessary  to  be  able  to  model  turbulent 
supercritical  flows.  Existing  models  such  as  the  k-e  turbulence  model 
may  require  additional  calibration  in  the  supercritical  regime.  Special 
considerations  must  be  given  to  the  validity  for  the  standard  wall 
functions.  Large  eddy  simulations  may  also  be  needed  for 
understanding  the  nature  of  turbulence  in  these  flows. 

d.  Analysis  of  Additives:  Anti-oxidants  are  added  to  the  fuel  to  inhibit 
thermal  oxidation  that  directly  contributes  to  the  formation  of 
insolubles.  Metal  deactivators  (MDAs)  have  also  been  used  to  reduce 
the  catalytic  activity  of  dissolved  metals.  However,  at  higher 
temperatures,  the  stability  of  these  additives  become  a  concern.  The 
chemistry  of  the  additives  has  to  be  imderstood  in  order  to  predict 
thermal  stability.  Accurate  models  have  to  be  developed  to  simulate 
the  effect  of  the  additives  over  a  range  of  conditions. 

e.  Validation  with  Benchmark  Data:  The  models  developed  during  this 
study  have  to  be  validated  against  benchmark  data.  Since  reliable  data 
for  many  fluids  in  the  supercritical  regime  are  almost  non-existent,  it  is 
necessary  to  design  an  experimental  setup  to  perform  appropriate 
experiments.  The  conditions  prevailing  inside  jet  fuel  systems  should 
be  reproduced  in  the  laboratory.  Experiments  need  to  be  carried  out  to 
measure  quantities  that  can  be  used  to  validate  the  computational 
models.  Work  is  underway  to  identify  candidate  experimental  sites  for 
the  Phase  n  work.  Some  of  these  are  (i)  University  of  Iowa,  (ii)  United 
Technologies  Research  Center  in  Hartford,  CT,  (iii)  Pratt  &  Whitney, 
FL,  and  (iv)  General  Electric  Aircraft  Engines  in  Cincinnatti,  OH. 

The  details  of  these  tasks  will  be  presented  in  the  Phase  n  proposal. 

7.4  Phase  III  Objectives 

For  CFDRC,  the  potential  applications  of  this  work  in  Phase  EQ  include: 
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a.  Commercialization  of  the  developed  software  for  analyzing  supercriti¬ 
cal  flows ; 

b.  Commercialization  of  customized  software  for  analysis  of  jet  fuel 
systems  with  capabilities  to  model  complex  geometries  such  as  heat 
exchangers,  injector  tubes,  coils,  etc.; 

c.  Analysis  and  development  of  design  correlations,  novel  concepts,  and 
devices  involving  applicatons  of  supercritical  flows;  and 

d.  Providing  technical  support  and  services  to  users. 

Effort  is  currently  underway  to  identify  would-be  commercial  users  of  this  product. 
Discussions  are  currenty  in  progress  with  Pratt  &  Whitney,  Delavan,  and  GE  Aircraft 
Engines  for  collaborative  work  during  Phase  n  and  for  potential  Phase  HI  (non- 
SBIR)  commitments. 
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